---
title             : "Income Inequality and Life Online: Natural Experiments on the Effects of Smartphone Payday Loan Ads on Psychological Stress"
author: 
  - name          : "Jihye Lee"
    affiliation   : "1 and 2"
    corresponding : yes   
    email         : "jihyelee@stanford.edu"
  - name          : "James T. Hamilton"
    affiliation   : "1"
  - name          : "Nilam Ram"
    affiliation   : "1 and 3"
  - name          : "Katherine Roehrick"
    affiliation   : "1"
  - name          : "Byron Reeves"
    affiliation   : "1"

affiliation:
  - id            : "1"
    institution   : "Department of Communication, Stanford University"
  - id            : "2"
    institution   : "Moody College of Communication, University of Texas at Austin"
  - id            : "3"
    institution   : "Department of Psychology, Stanford University"
output: html_document
---
This script conducts analyses for the manuscript entitled "Income Inequality and Life Online: Natural Experiments on the Effects of Smartphone Payday Loan Ads on Psychological Stress" at Information, Communication, & Society. Online first: https://www.tandfonline.com/doi/abs/10.1080/1369118X.2022.2109982?journalCode=rics20.

```{r, load_packages, include = FALSE}
library(dplyr)
library(ggplot2)
library(lubridate)
library(glmmTMB)
library(scales)
library(ggeffects)
library(performance)
library(data.table)
options(scipen=999)
'%notin%' <- Negate('%in%')
```

#Read the dataset and run multilevel analyses for each dependent variable.
```{r}
payday.df <- readRDS(file = "update your working directory path here/ICS_data_public.rds")

#######Future-oriented content
#First, test for zero inflation
check_zeroinflation(glmmTMB(focusfuture_raw_words ~ low.income * payday.exposure * social_safety_nets + age_in_2019 + sex + race_recoded1 + edu_recoded1 + (1 |subject/experiment), 
                            family = poisson(),
                            data = payday.df),
                            tolerance = 0.05)

#Function to transform variables for beta regression because values of 0 and/or 1 occur 
#See Cribari-Neto & Zeileis (2010)
transform_perc <- function(prob) {
    (prob * (length(prob) - 1) + 0.5)/length(prob)
}
payday.df$focusfuture_conv2 <- transform_perc(payday.df$focusfuture_conv)

#Run a multilevel model
outcome <- "focusfuture_conv2"
fixed_fx <- "~ low.income*payday.exposure*social_safety_nets + age_in_2019 + sex + race_recoded1 + edu_recoded1 "
random_fx <- " + (1|subject/experiment)"

focusfuture_beta_3way <- glmmTMB(as.formula(paste(outcome, fixed_fx, random_fx)),
                                family=beta_family(),
                                data = payday.df)
summary(focusfuture_beta_3way)

#Calculate conditional and marginal R-squared after correcting the random effect structure. 
random_fx2 <- " + (1|subject:experiment)"
focusfuture_beta_3way2 <- glmmTMB(as.formula(paste(outcome, fixed_fx, random_fx2)),
                                family=beta_family(),
                                data = payday.df)
r2(focusfuture_beta_3way2) 

#######App switching 
#First, test for zero inflation -> probable zero-inflation
check_zeroinflation(glmmTMB(my_adist ~ low.income * payday.exposure * social_safety_nets + age_in_2019 + sex + race_recoded1 + edu_recoded1 + (1 | subject/experiment), 
                            family = poisson(),
                            data = payday.df),
                            tolerance = 0.05)

# As we found excessive zeros in the data distribution, we used a zero-inflated beta generalized linear mixed model. 
payday.df$app_switch_conv2 <- payday.df$app_switch
payday.df$app_switch_conv2[payday.df$app_switch==1] <- 0.99

# Run a multilevel model 
outcome <- "app_switch_conv2"
fixed_fx <- "~ low.income*payday.exposure*social_safety_nets + age_in_2019 + sex + race_recoded1 + edu_recoded1"
random_fx <- " + (1|subject/experiment)"

app_zbeta_3way <- glmmTMB(as.formula(paste(outcome, fixed_fx, random_fx)),
                          ziformula = ~.,
                          family = beta_family(),
                          data = payday.df)
summary(app_zbeta_3way)

#Calculate conditional and marginal R-squared after correcting the random effect structure.
random_fx2 <- " + (1|subject:experiment)"
app_zbeta_3way2 <- glmmTMB(as.formula(paste(outcome, fixed_fx, random_fx2)),
                           ziformula = ~.,
                           family=beta_family(),
                          data = payday.df)
r2(app_zbeta_3way2) 

#######Negative emotional content
#First, test for zero inflation
check_zeroinflation(glmmTMB(negemo_raw_words ~ low.income * payday.exposure * social_safety_nets + age_in_2019 + sex + race_recoded1 + edu_recoded1 + (1 | subject/experiment), 
                            family = poisson(),
                            data = payday.df),
                            tolerance = 0.05)

#Transform variables for beta regression because values of 0 and/or 1 occur 
payday.df$negemo_conv2 <- transform_perc(payday.df$negemo_conv)

#Run a multilevel model
outcome <- "negemo_conv2"
fixed_fx <- "~ low.income*payday.exposure*social_safety_nets + age_in_2019 + sex + race_recoded1 + edu_recoded1"
random_fx <- " + (1|subject/experiment)"

negemo_beta_3way <- glmmTMB(as.formula(paste(outcome, fixed_fx, random_fx)),
                         family=beta_family(),
                         data = payday.df)
summary(negemo_beta_3way)

#Calculate conditional and marginal R-squared after correcting the random effect structure.
random_fx2 <- " + (1|subject:experiment)"
negemo_beta_3way2 <- glmmTMB(as.formula(paste(outcome, fixed_fx, random_fx2)),
                                family=beta_family(),
                                data = payday.df)
r2_nakagawa(negemo_beta_3way2)

#######Positive emotional content
#First, test for zero inflation
check_zeroinflation(glmmTMB(posemo_raw_words ~ low.income * payday.exposure * social_safety_nets + age_in_2019 + sex + race_recoded1 + edu_recoded1 + (1 | subject/experiment), 
                            family = poisson(),
                            data = payday.df),
                            tolerance = 0.05)

#Transform variables for beta regression because values of 0 and/or 1 occur 
payday.df$posemo_conv2 <- transform_perc(payday.df$posemo_conv)

#Run a multilevel model
outcome <- "posemo_conv2"
fixed_fx <- "~ low.income*payday.exposure*social_safety_nets + age_in_2019 + sex + race_recoded1 + edu_recoded1"
random_fx <- " + (1|subject/experiment)"

posemo_beta_3way <- glmmTMB(as.formula(paste(outcome, fixed_fx, random_fx)),
                         family=beta_family(),
                         data = payday.df)
summary(posemo_beta_3way)

#Calculate conditional and marginal R-squared after correcting the random effect structure.
random_fx2 <- " + (1|subject:experiment)"
posemo_beta_3way2 <- glmmTMB(as.formula(paste(outcome, fixed_fx, random_fx2)),
                                family=beta_family(),
                                data = payday.df)
r2_nakagawa(posemo_beta_3way2)

```

#Figure 3 in Main Text
```{r}
#######Predicted group differences in smartphone use during pre- and post- payday loan ad exposures: Future-oriented content
future_pr_3way <- ggpredict(focusfuture_beta_3way, c("payday.exposure", "low.income", "social_safety_nets"))  
focusfuture_plot_nosocial <- ggplot(future_pr_3way %>% filter(facet == "social safety net: no"), 
                              aes(x=x, y=100*predicted, fill = group, color = group)) + 
                              theme_classic() +
                            geom_point(position = position_dodge(width=0.1)) + 
                            geom_line(aes(x=as.numeric(x), y=100*predicted, colour = group), position = position_dodge(width=0.1)) +
                            geom_errorbar(aes(ymin=100*conf.low, ymax=100*conf.high, color = group), 
                                          position=position_dodge(0.1), 
                                          width=.1) +
                            scale_color_manual(values=c("black", "blue"))  +
                            scale_x_discrete(labels= c("Pre-Ad", "Post-Ad")) +
                            labs(title = "No Safety Net") +
                            ylab("Future-oriented words (%)") +
                            theme(plot.title = element_text(size=20, hjust = 0.5), 
                                                          axis.text=element_text(size=20),
                                                          axis.title=element_text(size=20), 
                                                          legend.text=element_text(size=20), 
                                                          strip.text.x = element_text(size = 20),
                                                          legend.title=element_text(size=20),
                                                         legend.position = "none",
                                                          axis.title.x = element_blank(),
                                                          panel.border = element_rect(colour = "black", fill=NA, size=1)) + 
                              scale_y_continuous(breaks = seq(0, 8, 2), limits = c(0, 8)) +
                              annotate("text", x=0.68, y=0.93, label= "Higher Income", size= 7.5) + 
                              annotate("text", x=0.74, y=2.5, label = "Lower Income", size= 7.5, color = "blue")
focusfuture_plot_nosocial

focusfuture_plot_social <- ggplot(future_pr_3way %>% filter(facet == "social safety net: yes"), 
                           aes(x=x, y=100*predicted, fill = group, color = group)) + 
                              theme_classic() +
                            geom_point(position = position_dodge(width=0.1)) + 
                            geom_line(aes(x=as.numeric(x), y=100*predicted, colour = group), position = position_dodge(width=0.1)) +
                            geom_errorbar(aes(ymin=100*conf.low, ymax=100*conf.high, color = group), 
                                          position=position_dodge(0.1), 
                                          width=.1) +
                            scale_color_manual(values=c("black", "blue"))  +
                            scale_x_discrete(labels= c("Pre-Ad", "Post-Ad")) +
                            labs(title = "Safety Net") +
                            ylab("Future-oriented words (%)") +
                            theme(plot.title = element_text(size=20, hjust = 0.5), 
                                                          axis.text=element_text(size=20),
                                                          axis.title=element_text(size=20), 
                                                          legend.text=element_text(size=20), 
                                                          strip.text.x = element_text(size = 20),
                                                          legend.title=element_text(size=20),
                                                         legend.position = "none",
                                                          axis.title.x = element_blank(),
                                                          panel.border = element_rect(colour = "black", fill=NA, size=1)) + 
                            scale_y_continuous(breaks = seq(0, 8, 2), limits = c(0, 8)) +
                            annotate("text", x=0.72, y=1.0, label= "Lower Income", size= 7.5, color = "blue") + 
                            annotate("text", x=0.68, y=2.2, label = "Higher Income", size= 7.5) 
focusfuture_plot_social

#######Predicted group differences in smartphone use during pre- and post- payday loan ad exposures: App switching
app_pr_3way <- ggpredict(app_zbeta_3way, c("payday.exposure", "low.income", "social_safety_nets"), type = "zero_inflated")
app_plot_nosocial <- ggplot(app_pr_3way %>% filter(facet == "social safety net: no"), 
                             aes(x=x, y=100*predicted, fill = group, color = group)) + 
                              theme_classic() +
                            geom_point(position = position_dodge(width=0.1)) + 
                            geom_line(aes(x=as.numeric(x), y=100*predicted, colour = group), position = position_dodge(width=0.1)) +
                            geom_errorbar(aes(ymin=100*conf.low, ymax=100*conf.high, color = group), 
                                          position=position_dodge(0.1), 
                                          width=.1) +
                            scale_color_manual(values=c("black", "blue"))  +
                            scale_x_discrete(labels= c("Pre-Ad", "Post-Ad")) +
                            labs(title = "No Safety Net") +
                            ylab("Application switching (%)") +
                            theme(plot.title = element_text(size=20, hjust = 0.5), 
                                                          axis.text=element_text(size=20),
                                                          axis.title=element_text(size=20), 
                                                          legend.text=element_text(size=20), 
                                                          strip.text.x = element_text(size = 20),
                                                          legend.title=element_text(size=20),
                                                         legend.position = "none",
                                                          axis.title.x = element_blank(),
                                                          panel.border = element_rect(colour = "black", fill=NA, size=1)) + 
                            scale_y_continuous(breaks = seq(0, 100, 20), limits = c(0, 100)) +
                            annotate("text", x=0.69, y=20, label= "Lower Income", size= 7.5, color = "blue") + 
                            annotate("text", x=0.75, y=5.2, label = "Higher Income", size= 7.5)
app_plot_nosocial

app_plot_social <- ggplot(app_pr_3way %>% filter(facet == "social safety net: yes"), 
                            aes(x=x, y=100*predicted, fill = group, color = group)) + 
                              theme_classic() +
                            geom_point(position = position_dodge(width=0.1)) + 
                            geom_line(aes(x=as.numeric(x), y=100*predicted, colour = group), position = position_dodge(width=0.1)) +
                            geom_errorbar(aes(ymin=100*conf.low, ymax=100*conf.high, color = group), 
                                          position=position_dodge(0.1), 
                                          width=.1) +
                            scale_color_manual(values=c("black", "blue"))  +
                            scale_x_discrete(labels= c("Pre-Ad", "Post-Ad")) +
                            labs(title = "Safety Net") +
                            ylab("Application switching (%)") +
                            theme(plot.title = element_text(size=20, hjust = 0.5), 
                                                          axis.text=element_text(size=20),
                                                          axis.title=element_text(size=20), 
                                                          legend.text=element_text(size=20), 
                                                          strip.text.x = element_text(size = 20),
                                                          legend.title=element_text(size=20),
                                                         legend.position = "none",
                                                          axis.title.x = element_blank(),
                                                          panel.border = element_rect(colour = "black", fill=NA, size=1)) + 
                            scale_y_continuous(breaks = seq(0, 100, 20), limits = c(0, 100)) +
                            annotate("text", x=0.69, y=15.5, label= "Higher Income", size= 7.5) + 
                            annotate("text", x=0.70, y=35, label = "Lower Income", size= 7.5, color = "blue")
app_plot_social

#######Predicted group differences in smartphone use during pre- and post- payday loan ad exposures: Negative emotional content 
negemo_pr_3way <- ggpredict(negemo_beta_3way, c("payday.exposure", "low.income", "social_safety_nets"))
negemo_plot_nosocial <- ggplot(negemo_pr_3way %>% filter(facet == "social safety net: no"), 
                            aes(x=x, y=100*predicted, fill = group, color = group)) + 
                              theme_classic() +
                            geom_point(position = position_dodge(width=0.1)) + 
                            geom_line(aes(x=as.numeric(x), y=100*predicted, colour = group), position = position_dodge(width=0.1)) +
                            geom_errorbar(aes(ymin=100*conf.low, ymax=100*conf.high, color = group), 
                                          position=position_dodge(0.1), 
                                          width=.1) +
                            scale_color_manual(values=c("black", "blue"))  +
                            scale_x_discrete(labels= c("Pre-Ad", "Post-Ad")) +
                            labs(title = "No Safety Net") +
                            ylab("Negative emotional words (%)") +
                            theme(plot.title = element_text(size=20, hjust = 0.5), 
                                                          axis.text=element_text(size=20),
                                                          axis.title=element_text(size=20), 
                                                          legend.text=element_text(size=20), 
                                                          strip.text.x = element_text(size = 20),
                                                          legend.title=element_text(size=20),
                                                         legend.position = "none",
                                                          axis.title.x = element_blank(),
                                                          panel.border = element_rect(colour = "black", fill=NA, size=1))+ 
                            scale_y_continuous(breaks = seq(0, 6, 2), limits = c(0, 6)) +
                            annotate("text", x=0.70, y=1.05, label= "Higher Income", size= 7.5) + 
                            annotate("text", x=0.75, y=1.95, label = "Lower Income", size= 7.5, color = "blue")
negemo_plot_nosocial

negemo_plot_social <- ggplot(negemo_pr_3way %>% filter(facet == "social safety net: yes"), 
                             aes(x=x, y=100*predicted, fill = group, color = group)) + 
                              theme_classic() +
                            geom_point(position = position_dodge(width=0.1)) + 
                            geom_line(aes(x=as.numeric(x), y=100*predicted, colour = group), position = position_dodge(width=0.1)) +
                            geom_errorbar(aes(ymin=100*conf.low, ymax=100*conf.high, color = group), 
                                          position=position_dodge(0.1), 
                                          width=.1) +
                            scale_color_manual(values=c("black", "blue"))  +
                            scale_x_discrete(labels= c("Pre-Ad", "Post-Ad")) +
                            labs(title = "Safety Net") +
                            ylab("Negative emotional words (%)") +
                            theme(plot.title = element_text(size=20, hjust = 0.5), 
                                                          axis.text=element_text(size=20),
                                                          axis.title=element_text(size=20), 
                                                          legend.text=element_text(size=20), 
                                                          strip.text.x = element_text(size = 20),
                                                          legend.title=element_text(size=20),
                                                         legend.position = "none",
                                                          axis.title.x = element_blank(),
                                                          panel.border = element_rect(colour = "black", fill=NA, size=1)) + 
                            scale_y_continuous(breaks = seq(0, 6, 2), limits = c(0, 6)) +
                            annotate("text", x=0.70, y=0.6, label= "Lower Income", size= 7.5, color = "blue") + 
                            annotate("text", x=0.73, y=0.9, label = "Higher Income", size= 7.5)
negemo_plot_social


#######Predicted group differences in smartphone use during pre- and post- payday loan ad exposures: Positive emotional content
posemo_pr_3way <- ggpredict(posemo_beta_3way, c("payday.exposure", "low.income", "social_safety_nets"))  
posemo_plot_nosocial <- ggplot(posemo_pr_3way %>% filter(facet == "social safety net: no"), 
                               aes(x=x, y=100*predicted, fill = group, color = group)) + 
                              theme_classic() +
                            geom_point(position = position_dodge(width=0.1)) + 
                            geom_line(aes(x=as.numeric(x), y=100*predicted, colour = group), position = position_dodge(width=0.1)) +
                            geom_errorbar(aes(ymin=100*conf.low, ymax=100*conf.high, color = group), 
                                          position=position_dodge(0.1), 
                                          width=.1) +
                            scale_color_manual(values=c("black", "blue"))  +
                            scale_x_discrete(labels= c("Pre-Ad", "Post-Ad")) +
                            labs(title = "No Safety Net") +
                            ylab("Positive emotional words (%)") +
                            theme(plot.title = element_text(size=20, hjust = 0.5), 
                                                          axis.text=element_text(size=20),
                                                          axis.title=element_text(size=20), 
                                                          legend.text=element_text(size=20), 
                                                          strip.text.x = element_text(size = 20),
                                                          legend.title=element_text(size=20),
                                                         legend.position = "none",
                                                          axis.title.x = element_blank(),
                                                          panel.border = element_rect(colour = "black", fill=NA, size=1))  + 
                              scale_y_continuous(breaks = seq(0, 20, 5), limits = c(0, 20)) +
                              annotate("text", x=0.70, y=5.7, label= "Higher Income", size= 6.5) + 
                              annotate("text", x=0.74, y=1.6, label = "Lower Income", size= 6.5, color = "blue")
posemo_plot_nosocial

posemo_plot_social <- ggplot(posemo_pr_3way %>% filter(facet == "social safety net: yes"), 
                           aes(x=x, y=100*predicted, fill = group, color = group)) + 
                              theme_classic() +
                            geom_point(position = position_dodge(width=0.1)) + 
                            geom_line(aes(x=as.numeric(x), y=100*predicted, colour = group), position = position_dodge(width=0.1)) +
                            geom_errorbar(aes(ymin=100*conf.low, ymax=100*conf.high, color = group), 
                                          position=position_dodge(0.1), 
                                          width=.1) +
                            scale_color_manual(values=c("black", "blue"))  +
                            scale_x_discrete(labels= c("Pre-Ad", "Post-Ad")) +
                            labs(title = "Safety Net") +
                            ylab("Positive emotional words (%)") +
                            theme(plot.title = element_text(size=20, hjust = 0.5), 
                                                          axis.text=element_text(size=20),
                                                          axis.title=element_text(size=20), 
                                                          legend.text=element_text(size=20), 
                                                          strip.text.x = element_text(size = 20),
                                                          legend.title=element_text(size=20),
                                                         legend.position = "none",
                                                          axis.title.x = element_blank(),
                                                          panel.border = element_rect(colour = "black", fill=NA, size=1)) + 
                            scale_y_continuous(breaks = seq(0, 20, 5), limits = c(0, 20)) +
                            annotate("text", x=0.74, y=2.2, label= "Lower Income", size= 6.5, color = "blue") + 
                            annotate("text", x=0.70, y=5.2, label = "Higher Income", size= 6.5)
posemo_plot_social
```